E-pile model of self-organized criticality 
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The concept of percolation is combined with a self-consistent treatment of the interaction between 
the dynamics on a lattice and the external drive. Such a treatment can provide a mechanism by which 
the system evolves to criticality without fine tuning, thus offering a route to self-organized criticality 
(SOC) which in many cases is more natural than the weak random drive combined with boundary 
loss/dissipation as used in standard sand-pile formulations. We introduce a new metaphor, the e- 
pile model, and a formalism for electric conduction in random media to compute critical exponents 
for such a system. Variations of the model apply to a number of other physical problems, such as 
electric plasma discharges, dielectric relaxation, and the dynamics of the Earth's magnetotail. 



A long-standing problem in the study of dissipative 
non-equilibrium dynamical systems is to explain how 
long-range temporal correlations with a "1//" power 
spectrum develop from local interactions as mass, en- 
tropy, and energy are exchanged with the environment. 
One approach to the problem can be found within a 
class of discrete lattice models which feature so-called 
self-organized criticality (SOC). The notion of SOC [l[ 
has emerged from the scale-free statistics generated from 
numerical simulations of the dynamics of such model sys- 
tems in response to a slow drive. A particularly clear ex- 
ample of SOC is the sand-pile model of which a few 
variants are known 0, H, 0, [l| . In sand-piles a small input 
perturbation triggers chain reactions of redistribution of 
a certain quantity (like mass or energy) as dictated by a 
prescribed thresholding condition such as a limit on occu- 
pancy per site. In response to the slow drive the system 
self-adjusts and evolves into a "critical" state without fine 
tuning of any external or control parameters. This idea 
of natural evolution to criticality via a self-organization 
process is central to all SOC models. The SOC state 
features no intrinsic length or time scale, its fluctuations 
are scale-free and are characterized by an inverse power- 
law power spectral density (PSD). By analogy with the 
traditional critical phenomena it has been argued that 
there is universal behavior, which could be described by 
a set of critical exponents and scaling relations. Crucial 
to SOC are the microscopic dynamical rules which deter- 
mine how the perturbations propagate across the system, 
and which define the universality class of the model Q . 
Several important ingredients and physical implications 
of SOC dynamics can be found in e.g., Refs. 0, H, @|. 

In this Letter we propose a self-consistent approach to 
SOC dynamics by including a feedback relation of the 
system at criticality on the drive. Such an approach 
provides a clear mechanism by which the system auto- 
matically evolves to criticality without fine tuning any 
of its parameters, and offers a route to SOC which in 
many physical settings is more natural than the weak ran- 
dom drive combined with the boundary loss/dissipation 
as used in standard sand-pile formulations. We analyze 



Q 







r 




















*• 


4 












M 










* 


* 


























1 












^ 










^ 












t 


► 










* 


* 










1 


► 






► 


























*- 







A(J) 



Storage 



c 



p =p<M) 



Release 



R 



Loading 



R 



ext 





r 








>* 






i 
































< 


y 












H 












M 




AT 






















Tl. 








► 




M 


M 


► 
























































.a 
























FIG. 1: The e-pile model: a dielectric formulation of the SOC 
properties and of the loading-storage-release cycle, generic of 
the slowly driven, dissipative, thresholded dynamical systems. 
Bottom: the changing condition of the lattice for the increas- 
ing bias, A<f>. Filled circles show nodes in the conducting 
state. (This figure is in color only in the electronic version.) 



the microscopic properties of the SOC state as a trans- 
port problem of electric charge. This approach differs 
from the conventional SOC models, in which one deals 
with the transport problem of charge-neutral particles 
whose dynamics is not influenced by an external field. 
We cast this charge-transport problem into a percolation 
problem, then combine the concept of percolation with a 
self-consistent treatment of the interaction between the 
dynamics on the lattice and the external drive. The SOC 
literature is rich in metaphors (sand-pile models, forest- 
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fire models, slider-block models, etc.). This is by some 
interpreted as a lack of clear concepts and rigor. We 
believe, however, that these metaphors are useful tools 
when a cross-diciplinary field in science is under concep- 
tual development. In this tradition we will introduce the 
e-pile model ("e" stands for "electric") as a metaphor 
for a model SOC system which combines percolation and 
self-consistent interaction with the external drive. 

In the e-pile model the dynamical rules of transport 
are dictated by the dielectric properties of the medium, 
and thus make more direct contact with the basic laws 
of particle motion and the polarization and conductivity 
response of the system. The model enables one to obtain 
the critical exponents from known properties of percola- 
tion [lO | an d of transport of charge in random/disordered 
media ll|, [H, [ljl, 14, The critical exponents are 

expressible in terms of the percolation indices and are 
numerically very close to values reported in Refs. [3, S] , 
although those works are based on different assumptions 
of the nature of the critical state and of the details of the 
interactions. We take this conformity as a manifestation 
of the inherent universality of the SOC state. 

Description of the model. — The model explores the 
properties of the transition from an insulating to a con- 
ducting state in a self-adjusting, disordered dielectric 
medium. Prospective applications of the model pertain 
to a variety of physical problems, such as electric plasma 
discharges, dielectric relaxation, and the dynamics of 
the Earth's magnetotail [III], but also by mathematical 
analogy in problems beyond electromagnetics. The SOC 
regime is achieved by locating the dielectric between the 
plates of a capacitor and applying an electromotive force 
as shown in Fig. 1. It is assumed that the dielectric con- 
sists of a homogenous d-dimensional lattice of N ~ L d 
nodes, where L is a characteristic system size (i.e., the 
distance between the plates). Each node can be in one 
of two states: (i) a ground state, in which it acts as an 
insulator; or (ii) an excited state, in which it acts as a 
conductor. The probability for a given node to be in 
the conducting state depends on the potential difference 
A<p across the capacitor and is given by a monotonically 
growing function p(A(f>) with p(0) = 0. Charge can only 
be conducted between nearest-neighbor nodes and only 
when these nodes are in the conducting state. 

The lattice as a whole acts as an insulator as long as 
there is no connected cluster of conducting nodes pro- 
viding at least one conduction path for electric charges 
across the lattice from one capacitor plate to the other. 
However, if the electromotive force V exceeds the critical 
value p(A(f> c ) — p c , where p c is the percolation thresh- 
old probability, the lattice spontaneously develops a con- 
ducting path and will as a whole suddenly become con- 
ducting, so that its resistance R drops from infinity to a 
finite value. As percolation is established, the capacitor 
starts to discharge, and we shall assume that the dis- 
charge time is short compared to the relaxation time of 
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FIG. 2: Dynamics of the critical state: cross-lattice electric 
current in a model in which the conducting network is fully 
reset after each discharge event. Adapted from Ref. [l8j]. 



the external circuit, i.e., CR <C Ci? cxt (see Fig. 1). When 
the potential difference A(f> drops below A(j> c the number 
of conducting nodes is reduced according to the proba- 
bility p(A<j>). The lattice then switches to the insulating 
state, and the cycle repeats. The final state of the system 
will be ever persisting charging-discharging fluctuations 
near the percolation point. These fluctuations are self- 
organized. The work presented in this Letter is based on 
the conjecture that these self-organized fluctuations are 
critical, i.e., scale- free. 

This conjecture can be verified by building a quantita- 
tive model which can be simulated on a computer. For 
this we need two additional ingredients: time will have 
to be discretized and we have to devise a model for how 
conducting nodes are randomly added or removed from 
the lattice as the probability p(Acj)) changes from one 
time step to the next. One conceivable model could be 
to select Np(A<fi) conducting nodes at random at every 
time-step. This, however, is hardly what we want from 
a SOC model, since the lattice will have no memory of 
its state at the previous time step. A more interesting 
model would be to add or remove NAp conducting nodes 
at random every time step. Here Ap is the change in oc- 
cupation probability from the running time step to the 
next. In this model the lattice remembers its present 
state when it moves to the next. 

Another ingredient we need is a model for the conduc- 
tivity across the lattice when it is in the conducting state. 
One approach to the problem is to assign a resistance 
and/or complex impedance to the connection between 
conducting nodes, and use Kirchhoff's equations. By ap- 
plying Kirchhoff's equations to the random resistor net- 
work one obtains its average frequency- and wave- vector- 
dependent conductivity, and then the dependence of the 
conductivity on the order parameter [171]. It is possible 
to realize a numerical model in which the total resistance 
over the lattice is calculated from Kirchhoff's rules and 
graph theory for the running distribution of conducting 
nodes. Results of this model in the Markovian (memo- 
ryless) limit when the entire lattice of conducting nodes 
is completely updated at random after each discharge 
event are summarized in Ref. 18|. In Fig. 2 we repro- 
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duce a plot from this study, showing the self-organized 
fluctuations of the electric current as the voltage drop 
Acj) fluctuates around the percolation threshold. 

Alternatively, one may consider to inject test-particles 
and allow them to hop between the nearest-neighbor con- 
ducting states in the direction of the voltage drop. This 
setting has important features in common with broadly 
studied directed-percolation and forest-fire models (Ref. 
Q for a brief overview). The analogy is particularly 
clear if one distinguishes three types of nodes: (i) empty 
nodes, i.e., nodes in the insulating or ground state; (ii) 
green trees, i.e., nodes in the conducting state, which for 
the time being do not contain particles; and (iii) burning 
trees, i.e., nodes in the conducting state which contain 
at least one particle. The charged particles are thought 
of as fire in this model. The difference is that the green 
trees are not burned by the fire, i.e., they remain green 
after having transmitted the charge, and their population 
is only affected by the environmental change in terms of 
the changing bias, A<j). There are variants of the model 
in which not more than one (fermion-type) or arbitrary 
many (boson-type) particles per node are allowed. 

The purpose of this Letter is, however, not to present 
a numerical study to prove the existence of critical be- 
havior, but rather to devise analytical methods to com- 
pute the critical exponents, assuming that the fluctua- 
tions are critical. Our analytical results conform with 
previously reported results of numerical simulation of 
sand-pile models 0, Q and with known estimates in the 
mean- field limit [8j. 

The percolation exponents. — Here for completeness 
of the discussion we briefly review the definitions of the 
percolation exponents v, (3, and p (Refs. [Io|, EH)- In 
the vicinity of the percolation point p c , the percolation 
correlation (the pair connectedness) length diverges as 
£ oc \p— p c \~ u for both p > p c and p < p c . For p > p c , the 
probability to belong to an infinite cluster is oc (p — 
p c )P and the dc conductivity of the cluster is a^ c oc (jp — 
Pc) 1 *- The infinite cluster is a fractal, with the Hausdorff 
dimension df — d — (3/ v. 

ac- conduction exponent. — If a fractal conducting net- 
work is exposed to an ac field, its frequency-dependent 
conductivity takes the form cr(u>) oc uf 1 where u) is the 
applied frequency. The ac-conduction exponent 77 has 
been obtained by Gefen et al. from a random- 

walk model of particle diffusion on percolation clusters: 
i] = /i/(2f + p, — f3). Using known estimates [lOj for v, 
f3, and p it is found that 77 ~ 0, 0.3, and 0.6 for d = 1, 2, 
and 3, respectively. The mean-field result, holding for 
d > 6, is 77 = 1. (Note that by mean-field we mean 
mean-field percolation, which is not to be confused with 
the mean-field theories of SOC Despite of some con- 
formity in the critical exponents there are subtleties of 
definition of the upper critical dimension: d u = 6 for 
percolation, and d u = 4 for SOC.) 

Power spectral density exponent. — This exponent is 



obtained from linear response theory. Let ~E(t, r) be the 
electric field at time t at point r in the bulk of the di- 
electric. The polarization response to this field is defined 
as P(t, r) = x(t ~ t')E(t', r)dt' where x(t - t') is a 
response function. In the frequency domain, P(w, r) = 
x(w)E(w, r) where x( w ) is the frequency-dependent com- 
plex susceptibility. The power spectral density (PSD) of 
the polarization field is given by S{u) = (|P(ct>, r)| 2 } = 
|x(^)| 2 (|E(w, r)| 2 ) where the angle brackets denote an 
ensemble average. One can see that |x(w)| 2 is the PSD 
of the polarization field when the driving electric field 
is an uncorrelated white noise signal. With use of the 
Kramers-Kronig relation x(w) oc P J cLo'o~(uj')/uj'(uj' — uj) 
it is found that x( w ) 0(1 w''^ 1 - Hence, S(ui) oc l/ui a with 
a = 2 - 2r). We have a ~ 2, 1.4, and 0.8 for d = 1, 2, 
and 3. The mean-field result is a = 0. Note that, in 1 
dimension the PSD of the SOC state conforms with the 
PSD of a Brownian random-walk process, i.e., a = 2. 
The result in 2 dimensions, i.e., a ~ 1.4, may be used 
to explain the low-frequency fluctuation spectrum in the 
current sheet of the Earth's magnetotail, which for years 
have been an issue in geo-space plasma research [16j . 

Dispersion-relation exponent. — The dispersion rela- 
tion shows how the lifetime of an "activation" cluster 
scales with its size [1| ■ In the e-pile model, by activation 
cluster one would mean a connected cluster of conducting 
nodes with particles. The finite lifetime of such clusters is 
due to diffusion of the particles in the conducting domain. 
We assume that the diffusion occurs as a result of hop- 
ping of charged particles between the nearest-neighbor 
conducting nodes and that there is a characteristic hop- 
ping time At after which a particle changes its site of 
residence by hopping at random from the hosting node 
to one of the nearest-neighbor conducting nodes. Given 
an activation cluster of linear size I one estimates its life- 
time to be t ~ iAt, where we have set the microscopic 
lattice distance to unity for simplicity. The hopping time 
is estimated as the microscopic conduction time, i.e., the 
inverse of the average size-dependent conductivity of the 
cluster, cr(£). The latter scales with the cluster size as 
a(£) oc [IH leading to t ~ £/a{l) oc Hence, 
the dispersion exponent is z = 1 + 77. It is found that 
z ~ 1.3 and 1.6 for d — 2 and 3, respectively. These 
values are numerically very close to the values obtained 
in the Zhang model [3j but with different analytical ex- 
pression. Similar values have also been found in the nu- 
merical simulation of Ref. . In the mean- field limit we 
have z = 2. This result conforms with the results of Ref. 
Q, in which a systematic mean-field treatment of SOC 
properties is given. In 1 dimension we have z = 1. 

Relaxation-function exponent. — The issue of the dis- 
persion relation is closely connected with the issue of re- 
laxation of charge-density inhomogeneities on a fractal 
lattice at percolation. Here we propose a self-consistent 
approach to the problem, in which the dynamics of decay 
of the activation clusters is governed by the electric field 
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fluctuations, produced by the electric charges themselves. 
Let pit, r) be the charge density at time t at point r in the 
bulk of the dielectric. These charges inside the dielectric 
generate the electric field inhomogeneity E(t, r) whose 
divergence is V • E(i, r) = 47r/3(t,r). The polarization 
response to this field is P(t, r) = x(t - *')E(*', r)dt' 
and its divergence is V • P(t,r) = —p(t,r). The flow 
function is defined by j(t, r) = d t P{t, r) where d t denotes 
time derivative. From the equation of continuity we find 
dtf>(t, r) = — 4-7T J 00 cr(t — t')p(t' , r)dt' where we have in- 
troduced a{t — t') = d t x(t — t')- If the dynamics start at 
time t = 0, then due to causality dt' = f Q dt' . One 
readily rewrites the dynamic equations in the Laplace do- 
main to find sp(s, r) + 4:ira(s)p(s, r) = p(0, r) [15j. Here 
p(0, r) is the initial charge density, which is set to 1 for 
simplicity. Assuming for the conductivity a(s) = ns 11 
with k a constant, then setting r\ = 1— £ we have p(s, r) = 
l/(s + Attks 1 ^^). The inverse Laplace transform of this 
is the Mittag-Leffler function [l9[ whose initial-time be- 
havior is a stretched-exponential decay function pit, r) oc 
exp [— 47TK^/r(C + 1)] . This stretched-exponential be- 
havior was suggested for SOC systems in Ref. 0], with- 
out a rigorous derivation. For the relaxation-function ex- 
ponent £ = 1 — rj we have £ ~ 0.7 for d = 2 and £ ~ 0.4 
for d = 3. The mean-field value is C = 0. Note that 
1-dimensional critical relaxation is exponential: £ = 1. 

It has been shown that the properties of ac conduc- 
tion and dielectric relaxation in disordered dielectrics are 
expressible in terms of fractional-derivative equations, 
and a systematic derivation of the fractional relaxation 
and fractional diffusion equations was presented in Ref. 

. These results suggest that systems in states of self- 
organized criticality are described by fractional kinetics. 

Size- distribution exponent. — In sand-pile models the 
exponent r is introduced to characterize the power-law 
distribution of avalanche sizes measured as the number 
of sites involved in the dynamics of isolated relaxation 
events 0, 0, S| ■ These notions are naturally present in e- 
piles, and the r values are readily obtained from Eq. (5) of 
Ref. 0] where one replaces the "noise" exponent <j> with 
a = 2(1 — rj), the fractal dimension D with df, and the 
dynamical exponent z with I + 77 to find (2/(3 — r) = a(l + 
77). The final result depends on the convention of order of 
limits. We have to specify which limit to be taken first: 
the thermodynamic limit L — > 00, or the zero-frequency 
limit ui — ► 0. If the thermodynamic limit is taken first, we 
make use of the percolation dimension df — d — [3/u and 
then apply r] = p/i2v+p—0) to obtain r ~ 2.0 for d = 2 
and r ~ 2.5 for d = 3. The mean-field result deriving 
from rj = 1 and a — is t = 3. If the zero-frequency 
limit is taken first, then one has to replace df with d, 
then to set 77 = in the dc (frequency-independent) limit 
to get z = 1, a — 2, and t = 3 — 2/d. This dependence 
of t on d was reported in Ref. 0]. Numerically, r ~ 2.0 
for d = 2 and t ~ 2.3 for d = 3. The corresponding 
mean- field value is obtained for df = 4 yielding r = 2.5, 



in agreement with Ref. [8 
t in Ref. [2| and in Refs 




Note that the definitions of 
differ by 1. 



Summary. — In conclusion, we have introduced a new 
self-consistent approach to realize SOC dynamics: the e- 
pile. We have provided analytical methods to calculate 
the main critical exponents characterizing the dynamics, 
under the assumption that the fluctuations are critical. 
Other critical exponents may be obtained if desired. The 
results agree with previously reported results from nu- 
merical investigations of the sand-pile models. Numer- 
ical simulation of the e-pile dynamics is under way for 
comparison with the analytical predictions. 
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